function [Xp, Yp, tOut, uOut] = Hydro2v3(Nx,Ny,Lx,Ly,tf,NSteps)
% Fourier Basis in Spatial Directions
%

global Step;

Step =0;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%% Initializations and Parameters                  %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Nt = Nx*Ny; % Total Number of sites in lattice
Nv = Nt*5; % Total number of variables

PlotData = 0;
delay = 1;

tSamples = 0:(tf/(NSteps-1)):tf;


% % % dt = 1.e-6;
% % % aTol = 1.e-7;
% % % rTol = 1.e-3;
% % % Fudge = 0.9;




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% **
    %%%%%%%%%    Build the Coordinate and Physical Grids      %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% % % [xx,XX,D1,D2]=BuildGrid('fourier','fourier');
% % % 
% % % x   = xx{1}; 
% % % X   = XX{1};
% % % %Dx  = D1{1};
% % % %Dxx = D2{1,1};
% % % %
% % % y   = xx{2}; 
% % % Y   = XX{2};
% % % %Dy  = D1{2};
% % % %Dyy = D2{2,2};
% % % %
% % % %Dxy = D2{1,2};
% % % 


Jx = 2*pi/Lx;
Jy = 2*pi/Ly;

x = -pi + 2*pi*(0:Nx-1)'/Nx;
y = -pi + 2*pi*(0:Ny-1)'/Ny;
[X,Y]=ndgrid(x,y);

xp = (Lx/(2*pi))*x;
yp = (Ly/(2*pi))*y;
[Xp,Yp]=ndgrid(xp,yp);

Grid = {Nx, Ny, Lx, Ly};


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%         Define Helper Macro Variables           %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%                 Reserve Storeage                %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% RK variables
u    = zeros(Nv,1);
uxy  = 0*X;
uyx  = 0*X;


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%      Set Variables to Initial Conditions        %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % Set all variables to zero, padded into sufficient space


t = 0;
iteration = 1;

To = 1; % T here is (4*pi/3)Temp
ampo = 0.01;
Tmax = To*(1+2*ampo);
Tmin = To*(1-2*ampo);

ux = 0*X;
uy = 0*X;
T = 0*X +To;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Random Noise in a fixed band
%
% for m = 3:9
%     phase = 2*pi*rand(1);
%     amp = ampo*(rand(1)-0.5);
%     T = T + amp*cos(m*Y + phase*cos(2*X));
%     phase = 2*pi*rand(1);
%     amp = ampo*(rand(1)-0.5);
%     T = T + amp*cos(m*X + phase*sin(2*Y));
% end
%

NumModes = 5;
ModeMin  = 60;
ModeBand = 4;
Amp0     = 0.1;
kDrive   = ModeMin+randi(ModeBand,NumModes,2);
aDrive   = Amp0 * rand(NumModes,1);
phDrive  = 2*pi*rand(NumModes,1);
phDrivet = 2*pi*rand(NumModes,1);
wDrive   = (1/8*pi)*rand(NumModes,1);

DriveFlag = 1;
DriveFunction = @(t) BandOfNoise(t,X,Y,NumModes, aDrive, kDrive, phDrive, phDrivet, wDrive);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Blowing in the wind

% To get behavior like louis, we need to take T idiculously small
T  = 0*X +1.0/40;
% Let's add some noise and cross-coupling, too...
ux = 0*X;% + 0.000001*sin(3*X+rand(1));
uy = 0*X + 0.00001*sin(7*X+rand(1));
P  = 0*X;
B  = 0*X;

u(    +1:1*Nt) = ux(:);
u(1*Nt+1:2*Nt) = uy(:);
u(2*Nt+1:3*Nt) = T(:);
u(3*Nt+1:4*Nt) = P(:);
u(4*Nt+1:5*Nt) = B(:);

DriveFlag = 1;
DriveFunction = @(t) DriveTwoJets(t,X,Y);

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Louis's Initial Conditions
%
% --- so in geometrized unitz computational_box = [0,10]x[0,10]
% --- initial density =1   (and the density goes like T^3)
% --- initial velocity profile -- > v_background^i = 0  (ie no boosted brane)
%     + perturbation^i = 0.5 * sin( l_p * Pi/10 * y)  \delta^i_x
% 
% with l_p, i've got a couple of cases =10 or =20  qualitatively they
% are the same, but things develop
% faster for larger l_p
% 
% I've evolved for T = 60  (and have others for much longer, higher
% resoln to see the asympt soln),
% the 'turbulent' behavior though takes less than this to show.%
%

% Louis's reported IC's 
T  = 0*X +1.0;
ux = 0*X + 0.5 * sin(5*Y);
uy = 0*X;
P  = 0*X;
B  = 0*X;

% Unfortunately, this doesn't match our results at all -- the whole fluid
% dissipates immediately long before anything turbulent obtains
% 
% To get behavior like louis, we need to take T idiculously small
T  = 0*X +1.0/50;
% Let's add some noise and cross-coupling, too...
ux = 0*X + 0.5 * sin(5*Y);
uy = 0*X + 0.0001*sin(3*X+rand(1)*cos(7*Y));
P  = 0.000001*sin(3*X+rand(1));
B  = 0.000001*sin(2*Y+rand(1));

u(    +1:1*Nt) = ux(:);
u(1*Nt+1:2*Nt) = uy(:);
u(2*Nt+1:3*Nt) = T(:);
u(3*Nt+1:4*Nt) = P(:);
u(4*Nt+1:5*Nt) = B(:);

DriveFlag = 0;

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%     Time-Evolve via Chebyshev/FFT Derivatives   %%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;

params = {Nx, Ny, (2*pi/Lx), (2*pi/Ly), DriveFlag, DriveFunction};

options = odeset('RelTol',1e-5,'AbsTol',1e-8,'OutputFcn',@(t,u,flag) RealTime(t,u,flag,params));

[tOut,uOut]=ode45(@(t,u) H2Derivs(t,u,params),tSamples,u,options);
%sol=ode15s(@(t,u) H2Derivs(t,u,params),[0 tf],u); SLOW!!!

toc


% % % % figure(1);
% % % % for n = 1:1:NSteps 
% % % %     
% % % %     uxo = reshape(uOut(n,      1:  Nx*Ny),Nx,Ny);
% % % %     uyo = reshape(uOut(n,Nx*Ny+1:2*Nx*Ny),Nx,Ny);
% % % %     
% % % %     uxy = Jy*SpectralGrad(uxo,'fourier',2);
% % % %     uyx = Jx*SpectralGrad(uyo,'fourier',1);    
% % % %     
% % % %     w = uxy - uyx;
% % % %     
% % % %     surf(Xp,Yp,w), shading interp, zlim([-0.4,0.4]), view(30,60) , light, lighting phong, shg, pause(0.05)
% % % % 
% % % %     n/NSteps
% % % %     
% % % % end


% % % % % Check constituitive relations
% % % % 
% % % % t=tOut(NSteps);
% % % % u=uOut(NSteps,:);
% % % % 
% % % % udot = H2Derivs(t,u,params);
% % % % 
% % % % ux = reshape(u(     1:  Nt),Nx,Ny); 
% % % % uy = reshape(u(1*Nt+1:2*Nt),Nx,Ny); 
% % % % T  = reshape(u(2*Nt+1:3*Nt),Nx,Ny); 
% % % % P  = reshape(u(3*Nt+1:4*Nt),Nx,Ny); 
% % % % B  = reshape(u(4*Nt+1:5*Nt),Nx,Ny); 
% % % % 
% % % % uxx = Jx*SpectralGrad(ux,'fourier',1);
% % % % uyx = Jx*SpectralGrad(uy,'fourier',1);
% % % % Tx  = Jx*SpectralGrad( T,'fourier',1);
% % % % Px  = Jx*SpectralGrad( P,'fourier',1);
% % % % Bx  = Jx*SpectralGrad( B,'fourier',1);
% % % % 
% % % % uxy = Jy*SpectralGrad(ux,'fourier',2);
% % % % uyy = Jy*SpectralGrad(uy,'fourier',2);
% % % % Ty  = Jy*SpectralGrad( T,'fourier',2);
% % % % Py  = Jy*SpectralGrad( P,'fourier',2);
% % % % By  = Jy*SpectralGrad( B,'fourier',2);
% % % % 
% % % % uxt = reshape(udot(     1:  Nt),Nx,Ny); 
% % % % uyt = reshape(udot(1*Nt+1:2*Nt),Nx,Ny); 
% % % % Tt  = reshape(udot(2*Nt+1:3*Nt),Nx,Ny); 
% % % % Pt  = reshape(udot(3*Nt+1:4*Nt),Nx,Ny); 
% % % % Bt  = reshape(udot(4*Nt+1:5*Nt),Nx,Ny); 
% % % % 
% % % % cc = ((sqrt(3)*pi) -(9*log(3)) +18)/18;
% % % % 
% % % % 
% % % % ConRel = P.*(2+ux.^2+uy.^2) + (1/2).*T.^(-2).*(1+ux.^2+uy.^2).^(-1).*(P.*((-2).*T.^2.*(2+ux.^4+ ...
% % % %   3.*uy.^2+uy.^4+ux.^2.*(3+2.*uy.^2))+4.*cc.*(Tx.*ux+Ty.*uy+Tt.*(1+ ...
% % % %   ux.^2+uy.^2).^(1/2)).*(2+ux.^4+3.*uy.^2+uy.^4+ux.^2.*(3+2.*uy.^2)) ...
% % % %   +(-1).*T.*((-2).*ux.*uy.*(ux.^2+(-1).*uy.^2).*((-1).*uxy.*(1+ ...
% % % %   uy.^2)+(-1).*uxt.*uy.*(1+ux.^2+uy.^2).^(1/2)+uyx+ux.^2.*uyx+ux.*(( ...
% % % %   -1).*uxx.*uy+(1+ux.^2+uy.^2).^(1/2).*uyt+uy.*uyy))+cc.*(uxx.*(2+ ...
% % % %   3.*uy.^2+uy.^4+ux.^4.*(3+2.*uy.^2)+ux.^2.*(3+(-2).*uy.^4))+(-1).* ...
% % % %   ux.*(2.*uxy.*uy.^5+uxt.*(1+ux.^2+uy.^2).^(1/2).*((-2)+uy.^2+2.* ...
% % % %   uy.^4))+2.*uy.*(1+ux.^2+uy.^2).^(1/2).*uyt+3.*uy.^3.*(1+ux.^2+ ...
% % % %   uy.^2).^(1/2).*uyt+(-2).*ux.^5.*uy.*uyx+ux.^3.*(2.*uxy.*uy.^3+ ...
% % % %   uxt.*(1+ux.^2+uy.^2).^(1/2).*(3+2.*uy.^2)+2.*uy.^3.*uyx)+2.*uyy+ ...
% % % %   3.*uy.^2.*uyy+3.*uy.^4.*uyy+ux.^4.*((-2).*uy.*(1+ux.^2+uy.^2).^( ...
% % % %   1/2).*uyt+uyy+(-2).*uy.^2.*uyy)+ux.^2.*((-1).*uy.*(1+ux.^2+uy.^2) ...
% % % %   .^(1/2).*uyt+2.*uy.^3.*(1+ux.^2+uy.^2).^(1/2).*uyt+3.*uyy+2.* ...
% % % %   uy.^4.*uyy))))+(-1).*T.*(T.^3.*(uxt.*uy.*(1+uy.^2).*(1+ux.^2+ ... 
% % % %   uy.^2).^(1/2)+uxy.*(1+uy.^2).*(1+ux.^2+uy.^2)+(1+ux.^2).*(ux.*(1+ ...
% % % %   ux.^2+uy.^2).^(1/2).*uyt+uyx+ux.^2.*uyx+uy.^2.*uyx))+4.*B.*(1+ ...
% % % %   ux.^2).*(1+uy.^2).*((-1).*uxy.*(1+uy.^2)+(-1).*uxt.*uy.*(1+ux.^2+ ...
% % % %   uy.^2).^(1/2)+uyx+ux.^2.*uyx+ux.*((-1).*uxx.*uy+(1+ux.^2+uy.^2).^( ...
% % % %   1/2).*uyt+uy.*uyy))+2.*cc.*(Px.*ux.*(2+ux.^4+3.*uy.^2+uy.^4+ ...
% % % %   ux.^2.*(3+2.*uy.^2))+(Py.*uy+Pt.*(1+ux.^2+uy.^2).^(1/2)).*(2+ ...
% % % %   ux.^4+3.*uy.^2+uy.^4+ux.^2.*(3+2.*uy.^2))+(-2).*B.*(uxy.*(1+uy.^2) ...
% % % %   .*((-1)+ux.^2.*((-1)+uy.^2))+(-1).*ux.^4.*uy.^2.*uyx+(1+uy.^2).* ...
% % % %   uyx+ux.^2.*(uxt.*uy.*(1+uy.^2).*(1+ux.^2+uy.^2).^(1/2)+uyx)+(-1).* ...
% % % %   ux.*uy.^2.*((1+ux.^2+uy.^2).^(1/2).*uyt+uy.*uyy)+ux.^3.*uy.*(uxx.* ...
% % % %   (1+uy.^2)+(-1).*uy.*((1+ux.^2+uy.^2).^(1/2).*uyt+uy.*uyy))))));
% % % % 
% % % % figure(1);
% % % % surf(Xp,Yp,ConRel), shading interp, view(30,60) , light, lighting phong, shg, pause(0.05)
% % % % 
% % % % abs(norm(ConRel(:)))/Nx*Ny

end










function status = RealTime(t,u,flag,params)

% % % length(t)
% % % t
global Step;

if strcmp(flag,'done')
    status = 0;
    return
end
    
Nx = params{1};
Ny = params{2};
Jx = params{3};
Jy = params{4};

Nt = Nx*Ny;

ux = reshape(u(     1:  Nt),Nx,Ny); 
uy = reshape(u(1*Nt+1:2*Nt),Nx,Ny); 

uyx = Jx*SpectralGrad(uy,'fourier',1);
uxy = Jy*SpectralGrad(ux,'fourier',2);

w = uxy - uyx;
wn = abs(norm(w(:)));


figure(1);
%            clf;
surf(w);
view(0,90),
%            zlim([Tmin Tmax])
title(['w(t) at t = ',num2str(t),'with norm(w) = ',num2str(wn)]);
shading interp;
%refresh();

Step=Step+1;
%saveas(gcf,['Hydro2v3_',num2str(Step)],'tif')

status = 0;
end


